home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-02 / polynom.zip / POLYNOM.PRG < prev   
Text File  |  1993-04-12  |  5KB  |  175 lines
  1. program PolyNom;  {Polynomial Curve Fitter pgm.}
  2.    uses crt;
  3.  
  4. const
  5.    RowSize   = 250;                                     {Mod. #1}
  6.    ColSize   = 2;
  7.    MaxDegree = 7;                                       {Mod. #2}
  8.  
  9. type
  10.    Array2Type     = array[1..RowSize, 1..ColSize] of real;
  11.    ArrayType      = array[1..RowSize] of real;
  12.    CoeffArrayType = array[0..30] of real;
  13.  
  14. var
  15.    DataArray               : Array2Type;
  16.    XArray, YArray          : ArrayType;
  17.    PowerArray, RHS, Coeffs : CoeffArrayType;
  18.    Factors                 : array[1..15, 1..15] of real;
  19.    NumDataPairs, Degree, TwoDegree, NumEqns, Code : integer;
  20.    J, K, L, M, Last                               : integer;
  21.    CharFlag                                       : char;
  22.    TimeToQuit, NoMoreToDo                         : boolean;
  23.    Work, Numer, Denom, Correlation, MeanY, X, Y   : real;
  24.  
  25. {$I GetNumI.PSL}
  26. {$I GetNumR.PSL}
  27. {$I Load2Arr.PSL}
  28. {$I Power.PSL}
  29. {$I Stats.PSL}
  30.  
  31. procedure Interpolate;
  32.  
  33. begin
  34.    NoMoreToDo := false;
  35.    repeat;
  36.       writeln;
  37.       writeln('X value to evaluate or E to End evaluations');
  38.       GetNumR(X, CharFlag, Code);
  39.       case Code of
  40.          -1: if CharFlag = 'E' then
  41.                 NoMoreToDo := true;
  42.           0: begin
  43.                 Y := 0.0;
  44.                 for K := 1 to NumEqns do
  45.                    Y := Y + Coeffs[K] * Power(X, K - 1);
  46.                 writeln('Y =', Y)
  47.              end
  48.       end
  49.    until
  50.       NoMoreToDo
  51. end;
  52.  
  53. procedure SolveEqns;
  54.  
  55. begin
  56.    TwoDegree := Degree * 2;
  57.    NumEqns   := Degree + 1;
  58.    for J := 0 to TwoDegree do
  59.       begin
  60.          PowerArray[J] := 0.0;
  61.          for K := 1 to NumDataPairs do
  62.             PowerArray[J] := PowerArray[J] + Power(XArray[K], J)
  63.       end;
  64.    for J := 1 to NumEqns do
  65.       begin
  66.          RHS[J] := 0.0;
  67.          for K := 1 to NumDataPairs do
  68.             RHS[J] := RHS[J] + YArray[K] *
  69.                                Power(XArray[K], J - 1)
  70.       end;
  71.    for J := 1 to NumEqns do
  72.       for K := 1 to NumEqns do
  73.          Factors[J,K] := PowerArray[J + K - 2];
  74.    for K := 1 to Degree do
  75.       begin
  76.          M := K + 1;
  77.          L := K;
  78.          repeat
  79.             if abs(Factors[M,K]) > abs(Factors[L,K]) then
  80.                L := M;
  81.             inc(M)
  82.          until
  83.             M > NumEqns;
  84.          for J := K to NumEqns do
  85.             begin
  86.                Work         := Factors[K,J];
  87.                Factors[K,J] := Factors[L,J];
  88.                Factors[L,J] := Work;
  89.             end;
  90.          Work := RHS[K];  RHS[K] := RHS[L];  RHS[L] := Work;
  91.          M := K + 1;
  92.          repeat
  93.             Work         := Factors[M,K] / Factors[K,K];
  94.             Factors[M,K] := 0.0;
  95.             for J := K + 1 to NumEqns do
  96.                Factors[M,J] := Factors[M,J] -
  97.                                Work * Factors[K,J];
  98.             RHS[M] := RHS[M] - Work * RHS[K];
  99.             inc(M)
  100.          until
  101.             M > NumEqns
  102.       end;
  103.    Coeffs[NumEqns] := RHS[NumEqns] / Factors[NumEqns,NumEqns];
  104.    for M := Degree downto 1 do
  105.       begin
  106.          Work := 0.0;
  107.          for J := M + 1 to NumEqns do
  108.             begin
  109.                Work      := Work + Factors[M,J] * Coeffs[J];
  110.                Coeffs[M] := (RHS[M] - Work) / Factors[M,M]
  111.             end
  112.       end;
  113.    writeln;
  114.    writeln('X Power', 'Coefficient':22);
  115.    for J := 1 to NumEqns do
  116.       writeln(J - 1:3, Coeffs[J]:28);
  117.    Numer := 0.0;
  118.    Denom := 0.0;
  119.    for J := 1 to NumDataPairs do
  120.       begin
  121.          Work := 0.0;
  122.          for K := 1 to NumEqns do
  123.             Work := Work + Coeffs[K] * Power(XArray[J], K - 1);
  124.          Numer := Numer + sqr(YArray[J] - Work);
  125.          Denom := Denom + sqr(YArray[J] - MeanY)
  126.       end;
  127.    if Denom = 0.0 then
  128.       Correlation := 1.0
  129.    else
  130.       Correlation := sqrt(1.0 - Numer / Denom);
  131.    writeln;
  132.    writeln('Correlation =', Correlation)
  133. end;
  134.  
  135. BEGIN
  136.    clrscr;
  137.    writeln('PolyFit - Polynomial curve fitter');
  138.    writeln;
  139.    writeln('The input data must now be read from a disk file.');
  140.    NumDataPairs := 0;
  141.    Load2Arr(DataArray, NumDataPairs, Code);
  142.    writeln;
  143.    if Code <> 0 then
  144.       begin
  145.          writeln('File was not read successfully.');
  146.          exit
  147.       end;
  148.    for J := 1 to NumDataPairs do
  149.       begin
  150.          XArray[J] := DataArray[J,1];
  151.          YArray[J] := DataArray[J,2]
  152.       end;
  153.    Stats(YArray, NumDataPairs, MeanY, Work, Work, Work, Work);
  154.    Last := NumDataPairs - 1;
  155.    if Last > MaxDegree then
  156.       Last := MaxDegree;
  157.    TimeToQuit := false;
  158.    repeat
  159.       writeln;
  160.       writeln('Degree to fit (0-', Last,
  161.                       ') or Q to Quit the program');
  162.       GetNumI(Degree, CharFlag, Code);
  163.       case Code of
  164.          -1: if CharFlag = 'Q' then
  165.                 TimeToQuit := true;
  166.           0: if (Degree >= 0) and (Degree <= Last) then
  167.                 begin
  168.                    SolveEqns;
  169.                    Interpolate
  170.                 end
  171.       end
  172.    until
  173.       TimeToQuit
  174. END.
  175.